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Abstract 

In this paper numerical simulations are used to calculate 
the turbulence dynamics simultaneously with the sound 
field for a high-speed near-sonic (Ma=0.9) compressible jet 
at two Reynolds numbers of 3,600 and 72,000. LES in 
conjunction with accurate numerical schemes is used to 
calculate the unsteady flow and sound in the near field of 
the jet. It is shown that the jet mean parameters, mean 
velocity fields and turbulence statistics arc in good 
agreement with experimental data and results from other 
simulations. The sound in the near-field is calculated 
directly from the simulations. The calculations are shown 
to capture the peak in the dilatation and pressure spectra 
around a Strouhal number St=0.25-0.3, in agreement with 
typical jet-noise spectra measured in experiments. 
Dilatation contours in the near-field show the formation of 
acoustic waves with a dominant wavelength of 3.2-4 jet 
diameters, corresponding to the peak in the dilatation 
spectra. As expected, the non-compact noise sources are 
found to be most dominant in the region corresponding to 
the end of the potential core. The contribution of the LES 
model to the radiated noise appears to be weak and does not 
contaminate the sound field with spurious high-frequency 
noise. However, the frequency spectra of the sound show a 
rapid fall-off away from the peak frequency. This is 
attributed to the quasi -laminar state of the shear- layers in 
the region prior to potential core closure, and a possible 
effect of insufficient azimuthal resolution at the observed 
location. Further analysis of the effect of the LES model, 
especially at high frequencies, is needed. 

1 Introduction 

The motivation of this work is to use LES techniques to 
calculate the acoustic emissions of jet engines. Once the 
jet simulations are validated in terms of turbulence 
dynamics, our focus is to predict the radiated jet noise, 
which is a dominant noise component for most aircraft jet 
engines at take-off conditions. A study of jet turbulence 
and its acoustics is relevant for many areas of applications 


including mixing enhancement for hot jet-exhaust plumes, 
but its primary motivation comes from the need to design 
more efficient engines with reduced noise emissions. 
Flight tests and wind-tunnel tests are useful but they are 
expensive, and quantitative measurements of those aspects 
of turbulence that represent the sources of noise radiation 
are very difficult, particularly in high-speed flows. This 
has led to the interest in using computational methods to 
try to better understand the noise generation. These 
insights should ultimately allow strategies for controlling 
or modifying the flow mechanisms to achieve the reduction 
of the noise emitted by jet engines. 

An accurate prediction of the sound radiated by a 
turbulent jet requires a method capable of reproducing the 
near-field turbulence dynamics with sufficient fidelity to 
allow the direct evaluation of the non-compact (distributed) 
sound sources. Obviously, RANS based methods require 
too much empirical input and are not suitable to describe 
accurately the distribution of the acoustic sources in space 
and time. A great deal of recent understanding of the 
turbulence physics in high-speed shear flows has come 
from direct numerical simulations (DNS) at low Reynolds 
numbers. In fact, Freund 11 has recently conducted DNS of a 
high-speed jet and its noise. However, as DNS is restricted 
to Reynolds numbers well below the values of engineering 
interest, LES techniques appear to be the only realistic 
available tool to obtain the necessary near-field flow data 
upon which to base the prediction of sound emitted by 
propulsive jets. Though still significantly more expensive 
than RANS methods, LES offers the advantage that little or 
no empirical input is needed, which is a significant 
advantage when one is interested in a robust method to 
predict the radiated sound field. This should allow us to 
better understand the role of the coherent structures to the 
noise generation. The- fact that the noise spectrum is 
dominated by the contribution of the large coherent 
structures justifies the use of LES for noise calculations. 
However, at very large Reynolds numbers the contribution 
of the smaller scales to the noise spectrum may be non- 
negligible in the range of the frequencies of interest, and 
this problem is yet to be investigated. 
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Over the last decade LES techniques have advanced to 
a point where they have been shown to predict complex 
flows (characterized by a large disparity between the 
different time and spatial turbulent scales) fairly accurately 
(e.g., Piomelli 2 and Wang 3 ). Use of LES for aeroacoustics 
applications is a natural development, but one which 
presents new challenges. The main challenge is related to 
the large disparity which exists between the energy levels 
associated with the fluctuations due to the large-scale 
turbulent motions and those associated with acoustic 
fluctuations. As most of the acoustic sources are situated 
in the early part of the turbulent region of the jet, a first 
requirement is to simulate correctly this turbulent flow so 
as to be able to calculate accurately the distribution and 
strength of the acoustic sources. Only the information 
about the high-frequency turbulence and hence its associated 
sound is lost in the LES method due to the inherent 
filtering. Furthermore, no simplifying assumptions about 
the acoustic wave propagation are required. The only 
difficulty is numerical; accurate simulation of sound waves 
concomitantly with the flow structures of the underlying 
turbulent flow is a considerable task. An important 
consequence is that accurate prediction of the radiated jet 
noise using LES techniques requires the use of numerical 
schemes with low dispersion and dissipation errors. 
Moreover, the quality of the noise data can be easily 
compromised by the boundary condition treatment. Higher 
Reynolds numbers and coarser grids put in turn a higher 
burden on the robustness and accuracy of the numerical 
methods used in aeroacoustic simulations. This is why 
attempts to simulate compressible jets and their radiated 
noise using LES are quite recent. 

Recently the group at NASA Lewis (e.g., Shih et 
al. 4 ) performed a three-dimensional LES calculation of a 
circular jet at Ma=L4 in a Ma=0.4 flight stream using a 
multiblock solver for curvilinear grids. A Kirchhoff 
method was used to propagate the acoustic emissions to the 
far-field. Choi et al. 5 simulated a Ma=1.2 jet with coflow 
at Ma=0.2, while Zhao et al. 6 have simulated a Ma=0.4 jet 
at Re=5,000. The success of most of these simulations 
was limited to a certain extent by the number of points 
used. This dictated what scales of fluctuations could be 
resolved and that in turn limited the frequency content of 
the simulated noise. Bogey et al. 7 conducted LES of a near 
sonic jet with a DRP scheme and artificial selective 
damping of high frequency waves on a Cartesian grid with 
6 million mesh points. They used the classical 
Smagorinsky model to calculate the subgrid (SGS) stresses. 
To decrease the cost of the computation the grid was refined 
only in a sector of the computational domain where data 
was collected in order to calculate the far-field radiated 
noise. Boersma and Lele 8 performed LES calculations of 
the near-field of a Ma=0.9 jet using a numerical method 
very similar to the one used here. The present work is a 
continuation of their efforts. Finally, it is relevant to 
mention efforts devoted to develop hybrid schemes which 
have a cost intermediate to LES and RANS. These 
schemes, though not as accurate as full LES computations 


of the near-field flow, are much less expensive. This 
approach is being developed at Penn State by Morris and 
co-workers (e.g., Morris et al. 9 * 10 ) and was applied to 
calculate circular and rectangular jets at Ma=2.I. In their 
method the Navier-Stokes equations were rewritten as 
equations for the nonlinear perturbations about the RANS 
solution to the mean flow. This made possible a 
decomposition of the instantaneous fluctuations into a 
time-averaged part, a resolved large-scale perturbation and 
an unresolved small-scale perturbation. The last part was 
accounted for using a SGS model. 

The paper is organized as follows. We start with a 
description of the numerical method, including the 
governing equations, implementation of the LES model, 
flow and boundary conditions, including the treatment of 
the flow equations at the polar axis. Next, we focus on the 
jet aerodynamics and we validate our simulations using 
experimental data obtained for turbulent jets at similar flow 
conditions. Finally, we use our simulations to investigate 
the radiated sound in the near-field of these jets. We 
conclude with a discussion of several issues that we aim to 
address in future work related to jet noise simulation using 
LES methods. 


2 Description of Numerical Method 
and Flow Conditions 

2.1 Numerical Method 

The primary focus of this paper is on near-sonic 
compressible cold jets with an acoustic Mach number 
Ma^U 0 /c 0 =0.9 and Reynolds numbers Re=U o (2r 0 )/v 
=3,600 and 72,000 based on the jet diameter at the inlet 
(2r 0 ), the jet centerline velocity at the inlet (U 0 ) and the 
speed of sound at infinity (c 0 ). At this Mach number the 
turbulent flow was investigated experimentally by Lau et 
al. 11 , and numerically by Freund 12 * 13 who performed a DNS 
simulation at Re=3,600, by Boersma and Lele 8 who 
investigated the flow at Re=3,600 and 36,000, and by 
Bogey et al. 7 who performed LES simulations at 
Re=65,000. These near-sonic conditions are motivated by 
the fan-stream exhaust conditions in a modem turbo-fan 
engine. Furthermore, due to the relatively high convection 
speed of the turbulent structures, the total computation 
time needed to obtain a statistically steady jet is also 
reduced. Our jet calculations are visualized in figure 1(a)- 
(b) using snapshots of the vorticity magnitude contours. 

The general numerical method is described in Freund 
et al. 14 , while Boersma and Lele 8 give the details of the 
implementation of the dynamic LES model in the original 
DNS code. Following Boersma and Lele 8 , the LES solver 
uses the non-density weighted compressible filtered 
variables (no Favre averaging). They found that the non- 
density weighted filtering as opposed to more widespread 
Favre filtering improves the robustness of the numerical 
solution, especially when compact schemes are used. They 
explained this by observing that the Favre averaged 
continuity equation is still a non-linear equation that can 
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Figure I: Instantaneous contours of vorticity magnitude, 
(a) Re=3,600 azimuthal section; (b) Re=72,000 azimuthal 
section; (c) Re=72,000 x=17.5r 0 ; (d) Re=72,000 x=3O.5r 0 


cause numerical instabilities with central-differences 
schemes. The downside to the non-density weighted 
filtering is that the form of the filtered equations is slightly 
more complex. The governing LES filtered equations are 
(an overbar is used to denote the grid filter): 

P.t +(pUi).i = “(puT “P=i>i (1) 

(pu,).. +(pu,u J ) J = 

_ _ __ ( 2 ) 

-P.i + T ,j.j -(pUiUj - pu.Uj) j - (pu; - pu,) , 

(E) , + ((E + p)^ ) } = -q, ■ + (TjjUj ) } 

__ (3) 

-(Euj - Euj + puj - pup j 


where Uj are the velocity components, p is the pressure, y 
is the ratio of specific heats, q ; =-kTj is the conductive 
heat-flux, f y = |i(Uj #j + u } s - 2/38^ k ) = 2pS£ is the 
Newtonian stress tensor, Sfj is the deviatoric part of the 
rate of strain tensor, E= p/(y - I) + O.SpUjUj is the resolved 
total energy density and p/(y-l) = pT/y is the equation 

of state. Note a further approximation related to using T 
rather than the Havre averaged temperature in the equation 
of state. The filtered continuity equation contains a subgrid 
mass flux while the Filtered momentum equation contains, 
besides the subgrid momentum flux, an additional term die 
to subgrid mass flux in the non-density weighted filtering. 
The presence of the subgrid mass flux in the continuity 
equation has the effect of improving the robustness 
properties of the numerical method. In equation (3) we 
neglected the unsteady term con taining the time variation of 
the subgrid kinetic energy (0.5(pUjUj -pu^)) as well as 
the additional subgrid term originating from filtering of 
TijUj. Using the equation of state and neglecting the 

convection of the subgrid kinetic energy by the resolved 
velocity, the last term in (3) repre senting the subgrid 
energy flux can be rewritten as (pUjT - pu-T). All these 

operators are discretized in cylindrical coordinates. The 
errors introduced by the non-commutativity of the filtering 
and discrete differentiation operations, as well as the 
contribution of the trace of the SGS stresses were also 
ignored. 

The code employs a centered six-order compact 
scheme to evaluate the spatial derivatives in the non- 
homogeneous directions, and Fourier spectral methods in 
the homogeneous (azimuthal) direction. The solution is 
advanced in time using a four-step Runge-Kutta method. 
These discretization schemes introduce very little artificial 
dissipation and allow sound waves to propagate accurately 
with only few grid points per wavelength (Freund and 
Lele 15 ). The number of modes is dropped near the polar axis 
(this is equivalent to Fourier filtering) so that the CFL 
constraint will be determined by the radial (or axial) 
spacing. This avoids the use of a very small time step in 
our explicit method. The time step is At^.Olro/Uo, 
corresponding approximately to CFL- 1. The computation 
was carried out on 32 processors of an 0rigin2000 using 
message passing interface (MPI). 

2.2 Subgrid Scale Model 

In LES the large, energy-containing scales are 
computed directly, while the small, unresolved scales that 
are nearly isotropic and their (non-linear) interaction with 
the large scales are modeled. However, there is no cut-off 
between the smaller scales and the larger ones. A lot of 
progress was made in LES of incompressible flows, 
particularly of jet flows. However, when applied to 
compressible flows especially for aeroacoustic applications 
the effect of the LES model on the flow fields is less 
understood. 
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In the present implementation, the LES filtering is 
implicitly defined by the computational grid used for the 
large-scale equations. The dynamic procedure is used to 
calculate the Smagorinsky-like constants in the SGS fluxes 
inequations (l)-(3). The main idea behind this approach is 
to determine dimensionless scaling coefficients in SGS 
models for the filtered nonlinear terms (Pierce 27 ). If such a 
term nfa*) is filtered (a/s are the dependent variables), it can 
be written as the sum of the resolved part and the modeled 
(subgrid) part: 

h(a i ) = n(a i ) + m(a i ) (4) 

In (4) the grid filter was denoted with an overbar and the 
filter width is A=^A Xf A^ A x ^ . A similar relation can be 

obtained if we apply a second filter, denoted by a hat 
symbol, and having a width A = 2A : 

ri(a i ) = n(a i ) + m(i i ) (5) 

Generally, the subgrid model is written as 
m(aj) = c s(a it A), where the ‘constant’ c now has a spatio- 
temporal variation. If we subtract the two previous 
relations, and allow the ‘constant’ c to pass unchanged 
through the test filtering operator, we obtain an equation in 
which aD terms are computable from the resolved field: 

n(ai ) - n(a ; ) = c (m(a t ) - m(a; )) (6) 

The left-hand side in (6) is called the Leonard term (L), 
while the right hand side is called the model term M. To 
obtain a single value for c (L and M are generally non- 
parallel vectors or tensors), equation (6) is renormalized and 
solved by least-squares following the procedure developed 
by Lilly 16 : 

c = (L'M)/(M’M) (7) 

where the square bracket denoted averaging in the 
homogeneous azimuthal direction. Only the ratio A/Ais 
needed in (7). 

Returning to the governing equation system (l M3), 
the subgrid terms in the continuity, momentum and energy 
equations are modeled according to (6). The dynamic 
coefficient corresponding to the subgrid mass, momentum 
and energy fluxes are calculated as follows: 

= , n j = p Uj , = A 2 |S|p j 

( 8 ) 

c s =(L ij M ij )/2(M ij M ij ) 

n ij = pUjUj ,m (j =pA 2 |S|S? j 


c E =-( L i M i )/(M 1 M i ) 

( 10 ) 

n i -pUjT , nij = pA 2 |S|Tj 

Note that, instead of calculating first the turbulent Prandtl 
number Pr T and then the equivalent ‘eddy-viscosity’ in the 
energy equation, v E = c E A J |s| = v T /Pr T , we calculated c E 
directly. In the previous relations ISI is the modulus of the 
strain-rate tensor Sj, and v T is the eddy viscosity in the 

momentum equations. The unsteady term in equation (2) 
does not require separate modeling; the subgrid mass flux 
model used in (1) is sufficient. 

To maintain stability of the solution field, the 
resolved flow variables are filtered in all three directions 
every 2vJq 0 time units (or every 200 time steps) using an 
explicit fourth-order accurate filter given by: 

= Mi-2 + 4f i-i + 10f| + f M - f i+2 )/ 16 (1 1) 

A new centerline treatment described in appendix A replaced 
the old method in which the equations were solved in 
Cartesian Coordinates. This allowed us to increase 
substantially (more than 12 times) the time interval at 
which this filter had to be applied. A detailed discussion of 
this issue is provided in Constantinescu and Lele 17 . As this 
filter is applied on the resolved variables, it is important to 
quantify its effect. In figure 2 the percentage change in the 
density field is shown. Note that the maximum change 
relative to the unfiltered density field is around 0.5%, with 
most of the change at much lower levels, including the area 
corresponding to the breakup of the potential core. Even 
though the effect of the damping on the large-scales 
structures is very low, there is still the possibility that this 
damping may affect the acoustic waves that originate in the 
near field. We expect that filtering the resolved variables 
can be avoided by using finer grids with near- unity aspect 
ratio ceils. In our simulation the aspect ratio is close to 
4.0 in most of the physical domain and this is believed to 
trigger some spurious oscillations. The development of 
higher-order schemes with behave robustly in treating 



Figure 2: Estimate of error introduced by filtering on the 
resolved density field 
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broadband non-linear problems (see, e.g. Nagarajan et al. 28 ) 
will eliminate the need for the introduction of any explicit 
numerical viscosity in the solution. 

The filter defined by ( 1 1) is also applied to the rate of 
strain tensor S t) at every time step before the terms in the 
SGS stresses are calculated This is done to avoid the 
formation of grid-to-grid oscillations in the solution due to 
under-resolved flow structures. This damping operation is 
important especially at the higher Reynolds number 
simulation where practically all the damping comes from 
the LES term. As noted by Boersma and Lele 8 this 
filtering affects only the sub-grid terms, and the error 
associated with it is smaller than the uncertainty involved 
with the use of the gradient hypothesis to obtain the SGS 
stresses. 

Snapshots of the SGS viscosity fields are displayed in 
figure 3(a)-{b), corresponding to the two Reynolds 
numbers. Again the difference in the position of the end of 
the potential core between the two calculations is obvious. 
However, if one takes this into account along with the 
scaling factor equal to the ratio between the two Reynolds 
numbers, the distribution of the SGS viscosity is very 
similar. Instantaneous values of the SGS viscosity relative 
to the molecular viscosity as high as 100 are observed at 
couple of points in the interior of the jet for Re=72,000, 
but ‘average’ values in these regions are about ten times 
less. 



Figure 3: Snapshots of eddy-viscosity fields normalized by 
molecular viscosity (a) Re=3,600; (b) Re=72,000 


2.3 Inflow Forcing 

As the shear layers are already turbulent at this 
Reynolds number, ideally we would need to have enough 
points to resolve the turbulence inside these shear layers. 
This would necessitate a lot of computational points inside 
the shear layers. The reader is referred to Freund and Lele 15 
for an in depth discussion of the issue. Additionally, the 
presence of the nozzle should be accounted in some way. 
All these requirements would increase substantially the 
computer requirements to perform such calculations. 
Instead, we started with laminar shear layers that were 
jittered to force transition. This is also evident from the 
snapshots of the vorticity magnitude shown in figure 1(a) 
and (b). The mean flow distribution at the inlet plane is 
assumed to be a rounded top-hat profile with periodic 
sinusoidal disturbances in the streamwise direction given 
by: 


u(r) = U 0 


f \ 1 J 

tanh 

2 2 


b 1 --*5 




n 


( 12 ) 


(I + a - sin(2re * St ♦ t )) 


where the Strouhal number St=2rof/U 0 is 0.9, the value of 
the thickness parameter is b=2.8 (corresponding to a 
momentum thickness of 0.09%) and the amplitude of the 
oscillations is a =0.005. Randomized azimuthal forcing 
with an amplitude of 0.025U o and zero mean is applied at 
the inlet plane to trigger the three-dimensional instabilities 
and finally transition to turbulence using the following 
distribution: 

u 9 =[o.025U 0 exp (-3(l-r/r 0 ) 2 )] 8 (13) 

where e is a random number between -0.5 and 0.5. The 
exponential function allows the random disturbances to be 
introduced only in the laminar shear layers. Low amplitude 
forcing was used because these disturbances do not satisfy 
the flow equations and generate spurious noise. This can 
be observed in the lower part of figure 9(a) where waves 
with wavelength close to 2r 0 are originating from the inlet 
region near the centerline. We could have forced earlier 
transition, especially at the lower Reynolds number, by 
increasing the amplitude of the disturbances, but we 
preferred to keep the spurious noise generated at the jet 
orifice to a small level. 


2.4 Boundary Conditions 

The computational domain extends to X rn =60% in the 
streamwise direction and R= 1 lr 0 in the radial direction. The 
computational grid consists of 320*192*64 points in the 
(x,r,9) directions, or a total of about 3.9 million mesh 
points, which is about an order of magnitude less than the 
one used by Freund 12 to calculate a similar jet at Re=3,600 
using DNS, but comparable to the meshes used by 
Boersma and Lele 8 and Bogey et a). 7 in their LES 
simulations. The points in the radial and streamwise 
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Figure 4: Variation of grid spacing in the computational 
domain (a) streamwise direction; (b) radial direction 


directions are clustered in the region r<5r 0 and 10r 0 <x<40r 0t 
where most of the sound sources are expected to be situated 
(see also figure 12(a)). In these regions the mesh spacing 
in the radial and axial directions is approximately 0.04r 0 and 
0. Ibr^ respectively. This mesh should be fine enough to 
allow an accurate representation of the large-scales 
structures. Figure 4 displays the distribution of Ax/r 0 and 
Ar/r 0 with the streamwise direction. One can see that near 
the boundaries the mesh spacing increases to 0.1 6r 0 and 
0.40^ respectively. The number of points in the radial 
direction corresponding to the initial core region (r<rg) is 
equal to 25, and the ratio between the initial momentum 
thickness and the radial grid spacing is 2.2. 

The frequency range of the fluctuations which can be 
captured on this mesh also needs to be kept in mind. Near 
r=0 the largest mesh spacing is Ax=0. 16r 0 . A wavelength 
with twice this spacing gives a Nyquist cut-off Strouhal 
number St c =6.9. However, as the sound waves spread 
outwards they encounter a mesh with increasing azimuthal 
spacing. At ^8,5^ a location where much of the acoustic 
data for the current LES is reported we have rA0=O.8r o 
which gives a cut-off Strouhal number St c =1.4. Evidently 
the range St>1.0 is poorly represented by the azimuthal 
mesh. This value of the Strouhal number is marked in the 
figures in which spectral information is presented at 
r=8.5r(). 

The formulation of the outflow and to a certain extent 
of the lateral boundary conditions for aeroacoustic 
computations is very important. This is because even a 
small reflection of a large-scale coherent structure exiting 
the computational domain can dominate the sound field, as 


generally the turbulence is several orders of magnitude more 
energetic than the sound. In the present simulations, zonal 
boundary conditions with artificial damping are used near 
the inlet and outlet to absorb outgoing disturbances before 
they reach the boundary and to avoid spurious noise 
generation via acoustic reflections at these boundaries. The 
target solution in these sponge layers was taken to be the 
self-similar solution of an incompressible jet. Non- 
reflecting boundary conditions are used at the lateral 
boundary, as well as damping terms in a layer close to this 
boundary to avoid the introduction of reflected waves into 
the domain. The inflow, outflow and lateral sponge layers 
are 2.5ro, 6.0^ and I.Or 0 wide, respectively. Taken 
together they contain roughly 20% of the computational 
points. The grid in the lateral and outflow sponge layers is 
stretched and an explicit 6 th order accurate filter is applied 
inside these layers every 40 time steps to further damp the 
near-grid scale turbulence and avoid reflections into the 
physical domain. 

3 Results 

3.1 Mean Properties and Turbulence 
Statistics of the Jet 

The general structure of the turbulent jet at is 
highlighted using snapshots of absolute vorticity contours 
in an azimuthal plane in Figures l(a)-Re=3,600 and (b)- 
Re=72,000 and in two plane situated at x/ro=17.5 and 30.0 
shown in figures 1(c) and (d) for Re=72,000. It is observed 
that vortical structures are generated in the initially laminar 
shear-layers due to the shear layer forcing. Quasi 
axi symmetric vortex rings are generated via Kelvin- 
Helmholtz instabilities as a result of the streamwise forcing 
of the m=0 mode at the inlet. Three-dimensionality leading 
to turbulence is accelerated by the interactions of these 
shear layers. The analysis of the jet mean velocity and 
turbulence statistics in Figures (5)-(6) shows that the 
transition towards turbulence starts at the end of the 
potential core situated around x=18r 0 for Re=3,600 and 
x=12r 0 for Re=72,000; the shear layers prior to the end of 
the ‘potential core* in these simulations are quasi- laminar. 
This has important consequences on the noise produced by 
this jet. The transition is followed by a region where the 
turbulence is fully developed and where the mean profiles 
will be self-similar in the streamwise direction. This 
region starts at x=24r 0 for the jet at Re=3,600, and x=18r 0 
for the jet at Re=72,000. 

We discuss in detail the distribution of the mean 
quantities and turbulence statistics for the jet at Re=3,6O0, 
but similar results were obtained for the other case, too. 
The region in which the jet is self-similar is found to 
extend up to x=40r 0 , after which the influence of the 
outflow sponge layer is felt. The statistics were collected 
over 8,000 time steps, corresponding to 80ro/c 0 time units. 
Before the statistics were collected, the simulation was run 
for about 18,000 time steps, corresponding to 3X n /c 0 ; all 
flow structures associated with the initialization of the 
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Figure 5: Streamwise evolution of: (a) Mean centerline 

velocity normalized by the inflow centerline velocity and 
its inverse; (b) Streamwise turbulent intensity; (c) the half- 
width of the jet normalized with the jet radius; Re=3*600 

simulation had enough time to be convected out of the 
domain. 

In figure 5(a)-(b) we show the distribution of the 
mean centerline velocity and the longitudinal turbulent 
intensities as a function of the distance to the jet orifice. 
Starting at the end of the potential core the centerline 
velocity decreases and some distance downstream the jet 
evolves to yield a self-similar jet. As the length of the 
potential core (core region where the centerline velocity is 
constant) is different in our simulation from the one 


corresponding to the experiment of Stromberg et a l. 18 as 
well as the DNS data of Freund 12 we shifted the data from 
these experiments/simulations by approximately 7^ 
corresponding to the difference in the length of the potential 
cores. The data of Lau et al. 11 was not shifted, even though 
the potential core in their experiments appears to end closer 
to llr 0 . The differences in the length of the potential cores 
and hence the position of virtual origin are not doe to a 
problem in our simulation* but rather they are a 
consequence of the absence of the jet nozzle in the 
simulations* and the shear-layer state (quasi-laminar). As 
already mentioned* we avoided the introduction of very 
strong disturbances in the flow that may have forced earlier 
transition* as these disturbances are known to radiate sound 
that will contaminate the sound radiated by the jet itself. 

Figure 5(a) clearly shows the decay of the centerline 
velocity with 1/x, as the inverse of the centerline velocity 
is seen to grow linearly with x* starting some short 
distance after the end of the potential core. One can also 
infer from the same plot the virtual origin of the jet 
corresponding to the intersection of that line with UqAJ^O* 
which is at x 0 =7.5r 0 . The decay constant 
B= U c (x)/U 0 ((x - x 0 )/2r 0 ) is reported to be between 5.4 

(Wygnanski and Fiedler 21 ) and 6.1 (Panchapakesan and 
Lumley 20 ), depending on the experimental setup. A least- 
square fit through our LES data gives a value for the decay 
constant B=5.9 which agrees well with the value found in 
the experiments of Hussein et al. 21 (B=5.9), the 
incompressible DNS simulation of Boersma et al. 23 
(B=5.9), and the LES simulation of Boerma and Lele 8 
(B=5.5). In fact, Boersma and Lele 8 suspected that the 
small length of their domain (45^ was responsible for the 
relatively lower value of the decay constant. Same cause 
may explain the results of Bogey et al. 7 who also found 
B=5.5 in their simulation. The length of their 
computational domain was 30^ compared with 6Or 0 in our 
simulations. 

The distribution of the streamwise non-dimensional 
velocity rms fluctuations a uu is presented in figure 5(b) 
along with the experimental data of Lau et ah 11 and 
Zaman 1 . The overall shape and range of turbulent 
intensities is correct* but the streamwise fluctuations are 
overestimated consistently compared to the two sets of 
experimental data* especially in the earlier part of the jet. 
The recent PIV data of Seiner et al. 29 is also plotted on the 
same figure for comparison. Unfortunately, it covers only 
the initial region of the jet core. We suspect some of the 
observed differences are due to the way in which transition 
to turbulence was forced in our simulations* via jittering* 
compared to their experiment where a nozzle was present. 
One should mention that Bogey 22 obtained even larger 
values for the longitudinal fluctuations in the fully 
developed region* while collecting statistics over a much 
longer time interval. Their data is also plotted for 
comparison in figure 5(b). Unfortunately, both Hussein et 
al. 21 and Panchapakesan and Lumley 20 do not provide this 
data. 
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Figure 6: Radial profiles of ihe Reynolds stresses in the 
self-similar region of the jet: (a) streamwise fluctuations 
O uu ; (b) radial fluctuations a vv ; (c) azimuthal fluctuations 
O ww ; (d) primary shear stress a uv ; R 0=3,600 


The growth of the jet is shown in figure 5(c) using 
the half-width of the jet 8 l/2 , defined as the distance at 

which the centerline velocity decays by a factor of two, 
plotted versus the streamwise coordinate. Again, we 
observe that starting with x/rp=24 in the fully-turbulent 
region the jet is growing linearly due to the lateral 
entrainment of potential flow from the irrotational region, 
and the slope of this line which defines the spreading rate of 
the jet is ^S, /2 /(x- x o )=0.096, a value that is close to 

the one deduced from the experiments of Panchapakesan et 
al. 20 (A=0.096) and Hussein et al. 21 (0.094) and the 
simulations of Boersma et al. 23 (A=0.095) and Bogey et al. 7 
(A=0.096). The location of the virtual origin deduced from 
the distribution of the half- width of the jet and the one 
inferred from the distribution of U</U c are very close, as 
expected. 

The radial profiles of the mean streamwise velocity, 
u, normalized by the centerline velocity, U c , are plotted 
versus the nondimensiona! radial coordinate B=r/(x-x 0 ) in 
figure 6(a) at five stations (x/r 0 =22, 26, 30, 34, 38) inside 
the self-similar region. As expected, these profiles are 
collapsing to a curve that is close to the gaussian velocity 
profile given by U(r)/U c = exp(-K(r/(x-x 0 )) 2 ), where 
K=ln2/A 2 =75.2. Panchapakesan and Lumley 20 also found 
K=75.2, while the simulations of Boersma et al. 23 and 
Bogey et al. 7 found K=76. 1 and K=75.2, respectively. 

The radial profiles of turbulent fluctuations u 2 , u 2 , 

Uq and primary shear stress u x u r in the self-similar region 
are also seen to collapse to curves that agree well with the 
experimental curve-fits (shown with symbols) obtained by 
Hussein et al. 21 and Panchapakesan and Lumley 20 . The 
non-dimensional turbulent intensities in the streamwise 
(a„ n ) and radial (a vv ) directions are shown in figure 6(b)-(d) 
versus T|, while the nondimensional shear stress (a uv ) is 
shown in Figure 6(e). The nondimensionaiization of the 
turbulent rms fluctuations and shear stresses is done using 

<J U . = ^|uj u j| /U c . The agreement appears to be better 

with the data of Hussein et al. for the streamwise turbulent 
intensities, while the radial intensities and the shear stress 
curves appear to be closer to the data of Panchapakesan and 
Lumley. The fact the agreement is not perfect may be 
caused by the insufficient radial extent of the computational 
domain, Reynolds number and Mach number differences 
between these experiments and our simulations, as well as 
the effect of the lateral and exit boundary conditions. 

Figure 7 shows power spectra of the resolved 
turbulent kinetic energy (normalized with the kinetic energy 
at the inflow) in the azimuthal direction at several 
streamwise locations corresponding to x/r^IO, 17, 24, 30, 
37 for r=r 0 r=2r ft and r=3ro, respectively. Except for stations 
at which the flow is not fully turbulent at that radial 
location, the form of the spectra shows that the flow is 
well resolved. It can be observed that the spectra at the 
stations situated far downstream (x/r 0 >24) have a slope that 
approaches -5/3 for the intermediate wavenumbers before 
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falling off at the high wavenumbers, in agreement with 
Kolmogorov’s theory. 



Figure 7: Azimuthal power spectra of the resolved kinetic 
energy at different stream wise locations for: (a) n=r 0 ; (b) 
r=2r 0 ; (c) r=4r 0 ; Re=3,600 


Finally, it is also of interest to look at the effect of 
the LES model on the jet turbulence. This is done by 
comparing at the relative strength of the different terms on 
the right-hand side of equation (2) for the streamwise 
velocity component for the Re=72,000 jet where we expect 
a larger effect of the LES model. In figure 8(a) we plotted 
the instantaneous fields corresponding to the sum of the 
convective, pressure and viscous terms, while in Figures 
8(b) and 8(c) we plotted the contributions from the sub-grid 


momentum flux and the time derivative of the subgrid mass 
flux, respectively. The contour levels go between 
±0.5p 0 c5/r 0 in figure 8(a), between ±O.01p o CQ/r 0 in 
figure 8(b), and between ±0.00 lp 0 CQ / r 0 in Figure 8(c). 
This shows that the LES contribution via the SGS stresses 
is far from being dominant. The results at Re=3,600 show 
a very similar picture in terms of relative magnitude of the 
LES terms relative to the resolved terms. In conclusion we 
can say that the agreement of the mean and turbulent 
properties between the present simulations and the available 
experimental data is satisfactory. It would be advantageous 
to compare this LES data with turbulence measurements in 
a Ma=0.9 jet. We are aware of several groups that are 
attempting such laboratory experiments, but published data 
are not available at the present time. 

3.2 Analysis of Near-Field Acoustic Data 

In this section we start using the LES fields as a 
database for acoustic calculations. A good measure for the 
sound waves emitted by a flow is the dilatation field. 
Snapshots of the dilatation fields in an azimuthal plane for 
both Reynolds numbers are shown in figures 9(a) and (b), 
while the fields in sections situated at streamwise locations 
defined by x/r 0 =23.5, 33, 36.5 and 42 are presented in 
figures 9(c)-(f), respectively. The presence of the patches of 
high positive and negative dilatation near the jet center that 
form a rather regular pattern corresponds to the large-scale 
vortex rings that are shed at the inlet excitation frequency 
(St=0.9). As these highly coherent structures start to 
interact, transition to turbulence take place and the 
distribution of the large scales becomes random. It is in 
this region that we expect most of the noise to originate. 

Away from the region where the noise sources are 
situated the dilatation is a direct measure of the pressure 
variation in time, in fact, the linearized (in vise id) energy 
equation gives: 

V-u = l —p t (14) 

Po c o 

We checked relation (14) at points situated away from the 
region where the non-linear interactions are expected to be 
important. For instance, the time series corresponding to 
the RHS and LHS in equation (14) nondimensionalized by 
r</c 0 are shown for a point situated at r=8.5r 0 and x= 25r 0 . 
The two time series are practically identical, when 
superimposed one on top of the other in figure 10(a). The 
pressure power spectrum at this point calculated starting 
from the pressure signal and the dilatation signal (obtained 
using the equivalent of (14) in Fourier space) are also very 
similar for the range of resolved frequencies as shown in 
figure 10(b). The differences are due to low -frequencies 
variations that are still present in the flow at this region, 
but it is clear that the dilatation oscillations at this location 
correspond to sound waves. 

With this in mind, one can try to describe 
qualitatively describe the sound generation using the 
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Figure 8: Residual in the streamwise momentum equationfor Re=72,000. Contribution from: 
(a) resolved terms Q RHS ; (b) sub-grid momentum flux Q^; (c) sub-grid mass flux 


dilatation fields. Figures 9(a) and (b) clearly show 
dilatation waves with a wavelength of about 6-8r 0 that are 
generated around a virtual origin situated at about x/^25 
for the Re=3,600 simulation, and at x/ro=20 for the 
Re=72,000 simulation. These locations correspond 
physically with the termination of the potential core, where 
most of the noise is expected to originate according to both 
experimental (Juve et al. 24 ) and DNS data (Freund 12 ). 

As in figure 9(a) the physical domain defined by 
(r>5r 0 , x>20r ft ) appears in to be in the acoustic ‘near-field* 
of the jet, meaning that the dilatation waves visible in the 
instantaneous dilatation field are traveling as sound waves, 
it is interesting to look at the time series of dilatation and 
pressure (figure 1 1) at a particular location corresponding to 
point P situated at r=8.5r 0 and x=36r 0 in figure 9(a). Indeed, 
these time series plotted in figure 1 1 clearly show the 
presence of a low-frequency component with a period of 
about 700 time steps (or 7.88^0^ corresponding to 
St-0.26. This Strouhal number matches the dominant 
high wavelength visible in the dilatation field in figure 
9(a). On top of the low-frequency waves, oscillations with 
a higher frequency in the range St= 1-4 are also observed. It 
is not clear if most of these higher frequency fluctuations 
are spurious. We suspect that these fluctuations can be 
reduced by running the simulation on a mesh with a 


smaller aspect ratio Ax/Ar (=4.0 in the present 
simulations). This is a task for future work. However, 
physically generated acoustic waves in this frequency range 
are expected to arise from turbulent shear layers. The 
presence of dominant low frequency waves at St=0.25-0.3 
also explains the dilatation patterns seen in figures 9(dH0, 
where these waves are traveling from the jet center toward 
the lateral boundary of the physical domain. Note that in 
these figures we represented the dilatation field over exactly 
this region that extends only lOr^ while the wavelength 
associated with the low frequency wave is around 7-8r 0 . 
This explain why in figure 9(d) and (f) we see only a region 
of positive dilatation in most of the domain, while the 
opposite is true in section 9(e) situated in between the 
previous two. These regions correspond to the 
compression and rarefaction fronts of the low-frequency 
wave that is seen in the dilatation time series in figure 11. 

In the lower part of figures 9(a) and 9(b) one can see 
the effect of the disturbances introduced at the jet inlet in 
terms of sound emissions. Though these disturbances 
radiate sound at the main excitation frequency (wavelength 
is close to 2r 0 ), they do not contaminate the sound radiated 
by the jet itself, which seems to be emitted at a much 
higher wavelength. The other important observation is that 
the outflow boundary appears to damp the turbulence 
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Figure 9: Instantaneous contours of dilatation magnitude, 4 contours equally spaced between +0.003c ( Vr Q: 

(a) Re=3,600 azimuthal section; (b) Re=72,000 azimuthal section; (c) Re=3,600 x=23.5r 0 ; (d) Re=3,600 x=33r 0 ; 
(e) Re=3,600 x=36.5r 0 ; (f) Re=3,600 x=42r 0 


without significant production of spurious sound that will 
propagate back into the domain and contaminate the sound 
pattern. 

We already alluded to the fact that that as the 
dominant sound waves originate from a region that 
corresponds to the termination of the potential core, the 
‘sound sources’ would be expected to be located in the same 
area. Further evidence for this is given by figure 12(a) 
where the instantaneous sound sources Srhs (except the 
contribution from the LES model denoted S LES , given in 
Figure 12(b)) are shown for the simulation at Re=72,000. 
The acoustic sources are calculated directly from their 
definitions: 


Srhs ~ c o(P u » u j ~ + (P c oP) u ^ ^ 

S L E S =Co(pu i u j -pu i u j ) i j 

Indeed, for this simulation the instantaneous acoustic 
sources of relatively high intensity are mostly located in 
the region 10r 0 <x<22r 0 and r<4r 0 . Their pattern and overall 
intensity levels (six equally spaced levels between 
±2.5pUo/ro) is quite close to that shown by Freund 12 . 
Comparison of figures 9(a) and 9(b) also suggests that the 
LES model does not become acoustically dominant, which 
is an important concern for application of LES to 
aeroacoustics. As the levels in figure 9(b) are between 
±0.25pUo /T q the dominant instantaneous acoustic sources 
are coming from the resolved -scale terms, while the 
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contribution of the LES term is at least one order of 
magnitude lower. 




Figure 10: Spectral content of the signal at a point 

situated at (p=8.5r„, x=25r 0 , 0=0) for the Re=3,600 jet (a) 
Time evolution of the dilatation and pressure time 
derivative, corresponding to the two terms in equation (14); 
(b) Pressure power spectra from the pressure and dilatation 
signals 
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Figure II: Time evolution of the pressure and dilatation at 
a point situated in the acoustic ‘near-field’ at (r=8.5r„, 
x=36.5r 0 , 0=0) for Re=3,600 jet 


Figures 13(a)-(c) show power spectra of the dilatation 
at points situated in the acoustic near field at r=8.5r 0 and 
x= 50rQ, 43r 0 and 33r 0 . These locations correspond to points 


SI, S2 and S3 in figure 9(a). Their distance from the 
‘virtual’ source center is 33^ 26r„ and 17r 0 for Re=3,600 
and at 39 r^, 32r^ and 23r^ for Re=72,000. The polar angles 
are 15°, 19° and 30° for Re=3.600 and 13", 15° and 22° for 
Re=72,000, respectively. The spectra are shown for both 
simulations, with the data for Re=72,000 being shifted 
down by a factor of 100 on the vertical axis to allow a 
better comparison. The noise spectra show that the range 
of Strouhal numbers corresponding to the peak noise 
emission is between 0.2 and 0.5, as expected. The power 
spectra at the higher Reynolds number appear to be more 
physical in the sense that the second peak at St=[-5 is 
much weaker and the shape of the spectra is closer to 
experimental results. Overall, one can observe the fact that 
for Re=3,600 the peak in these spectra is found around 
St=0.25 while the spectra at Re=72,000 appear to peak at a 
slightly higher value of St=0.3. This is evident if vre 
compare the spectra at points SI in the Re=3,600 jet aid 
S2 in the Re=72,000 jet that are situated approximately at 
the same radius from the end of the potential core and have 
same directivity angle. The second value (given in terms of 
dominant wavelength X=6.6r 0 ) agrees very well with the 
value found in the simulation of Bogey et al 7 at 
Re=68,000. 



L Re=72,000 
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Figure 12: Instantaneous visualization of the Lighthill 

sound source S using 6 equally spaced contours between 
±2.5pU 0 /r 0 . Positive levels are solid, negative levels are 
dotted, (a) contribution of the resolved variables S^; (b) 
LES sub-grid contribution S LES for Re=72,000 jet 
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The dilatation contours in the source region 
sometimes show the formation of a series of wave groups 
that seem to radiate sound very intensely for a certain period 
of time. These wave packets with a wavelength smaller 
than 2r 0 seem to form and be destroyed on a relatively short 
time-scale in a rather random fashion. They appear to be 
correlated with the relative peak in the noise spectra at 
St>l. However, as the cut-off Strouhal number at 
locations where the acoustic near-field is measures is close 
to these values, one probably should discard all together the 
spectral information for St> 1 .0. This is why a thick line 
is drawn on our spectra at St=1.0. 


icrV 



Finally, we use the pressure time series at points S I- 
S3 to calculate the sound pressure levels (SPL) in decibels 
for the simulation at Re=72,000. As observed in figure 14, 
the maximum sound pressure level is obtained around 150- 
!60dB. which is somewhat high but this is expected since 
the points are situated only 8.5r 0 from the jet axis. The 
overall sound pressure levels (OASPL) for points S 1 -S3 are 
129, 134 and 133dB, respectively. The decay in the SPL 
between St=0.2 and St= 1 is about 28dB (figure 14(b)), 
which is rather sharp as in translates into a decay of the 
SPL with frequency following a power law with n=-4, 
while experimental data suggest n=-2.5 to -2. However, 
one may speculate that the decay at higher frequencies is 
sharp because in our LES simulations we are not capturing 
the higher frequency sound produced by the turbulent shear 
layers prior to the potential core closure. Recall that these 
shear layers were quasi -laminar. There is also some effect 
of the increasing azimuthal spacing as St= 1.0 is 
approached. Better resolved near-field data is needed to sort 
out these effects. At the other end of the spectrum our 
results do not seem to show the very energetic low 
frequency waves that were present in the computed spectra 
of Bogey et al. 7 . This may be a consequence of the longer 
domain and somewhat different treatment of the equations 
near the outlet boundary that were suspected to be the cause 
of these waves. 




Figure 13: Dilatation power spectra for the simulations at 
Re =3, 600 and Re=72,000 at points (a) SI (^8.5^ x=50r 0 , 
9=0); (b) S2 (r=8.5ro, x=43r 0 , 0=0); (c) S3 (r=8.5r„, x=33r 0 
, 9 = 0 ) 


Figure 14: Sound pressure level obtained from the pressure 
at the points SI, S2 and S3 defined in figure 13 for the 
simulation at Re=72,000 (a) Strouhal number range 
between 0.01 and 5.0; (b) detail of the Strouhal number 
range between 0.01 and 1.5 


13 

American Institute of Aeronautics and Astronautics 






4 Conclusions 

In this paper we described the results of numerical 
simulations of two jets at same Mach number Ma=0.9 but 
at different Reynolds number of Re=3,600 and Re=72,000 
with the aim of showing the feasibility of calculating 
directly the sound sources and near-field noise using LES. 
The jet was excited randomly at the inlet plane to force 
transition. The code uses sixth-order compact schemes to 
evaluate the derivatives in the radial and stream wise 
directions while the evaluation of the derivatives in the 
azimuthal direction is done in the Fourier space. This 
ensures that a very little amount of artificial dissipation is 
introduced, allowing the evaluation of the influence of the 
SGS model separately from the effect of numerical 
dissipation, as well as making the code suitable to compute 
directly the radiated noise along with its aerodynamic 
properties. 

The First part of the paper dealt with the description of 
the jet aerodynamic flow characteristics in terms of the 
main jet parameters that characterize the growth in the fully 
developed region, mean profiles, and turbulence statistics. 
The simulation in terms of aerodynamic data were validated 
successfully by comparing with available experimental data 
as well as results of a DNS simulation carried out by 
Freund 12 at Re=3,600, as well a recent LES simulation 
carried out by Bogey et ai. 7 at Re=68,000. Next the noise 
computed directly in the near-field region corresponding to 
the physical domain of our simulation was investigated. 
As expected, the sound sources were found to be situated in 
the region near the end of the potential core, and the 
formation of sound waves was captured. The power spectra 
of these sound waves had its peak around a Strouhal 
number of 0.25-0.3, in agreement with various 
experimental studies. A more extensive validation of the 
sound results is under way. All these results establish that 
we have a numerical method that allows us to investigate 
in details the different mechanisms of sound generation. 

An ongoing goal of this work is to establish a 
benchmark LES database for cold jets, that would contain 
similar calculations at higher Reynolds number (in the 
range of 10 6 ), with the final aim to determine the range of 
frequencies over which reliable noise data can be extracted 
from LES and the range of frequencies that would require 
modeling. We intend to use the near-field data to compute 
the far-fiekl sound using either a Kirchoff integral method 
or a wave equation solver. The further development of 
hybrid methods to calculate the far-field radiated noise 
starting from direct calculations of the sound in the near- 
field is essential, especially as more realistic exhaust jet 
engine configurations are considered. 

Our preliminary results suggest that the LES model 
contribution to the radiated sound is not significant. To 
address this issue in more details we intend to use the DNS 
database of Freund 13 to extract space-time correlations of 
acoustic sources and compare with our LES results) at the 
same Reynolds number. This would allow to 
quantitatively estimate the attenuation and suppression of 
sound, especially at high frequency due to the unrepresented 


scales in LES, as well as to investigate in more details the 
noise contribution of the LES model. We estimate that a 
more consistent analysis on how the quality of jet-noise 
predictions depends on the mesh sizes, filtering, etc is 
necessary to further demonstrate the feasibility of LES 
methods to aeroacoustics applications. 

The question of estimating the higher-frequency sound 
information that is missing from our simulations because 
on one hand the LES filtered equations were used, and on 
the other due to quasi-laminar shear layers, requires further 
investigation. As pointed out by Freund and Lele 15 
predictions methods that allow the acoustic output of the 
small scales to be estimated and combined with the acoustic 
output of the resolved scales in LES are highly desirable. 
We expect that this database will provide valuable help in 
efforts related to the development of SGS acoustic models 
that are necessary to account for the noise information lost 
due to the LES filtering and, in general, the development of 
LES techniques tailored for noise predictions. 
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Appendix A: Centerline treatment 

As the governing system of equations is discretized in 
cylindrical coordinates, special care is given to the 
treatment of the equations at the polar axis, due to the 
presence of singular terms at r=0. Our experience shows 
that the quality of the LES solution, especially when 
compact finite-differences schemes are used, is especially 
sensitive to the type of equation treatment at the singularity 
axis. 

In the present work we use a new formulation in 
which a set of exact equations at the singularity axis is 
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derived using the appropriate series expansions for the 
variables in the original set of equations. The main idea is 
to reinterpret the regularity conditions developed in the 
context of pseudo- spectral methods. Besides increasing 
considerably the robustness of the numerical method 
compared to previous versions of the code, an advantage of 
the proposed treatment is that it preserves the same level of 
accuracy as for the interior scheme. 

The reader is referred to Constantinescu and Lele 17 for 
a detailed description of the method. Here, we will 
emphasize only the main points. The governing system of 
equations ( l)-(3) can be written in compact form as: 


^=RHS(Q) 

at 


(Al) 


where in our case the vector of unknowns is 
(Xpu^pi^pu^p.e) an d the right-hand-side term (RHS) 
contains the usual operators in cylindrical coordinates 
associated with the continuity (1), momentum (2) and 
energy equations (3), including the LES terms. Following 
Boyd 25 , the most general expansion of a single-valued " 
quantity (S) at the polar axis can be written as: 


S(r, 0) = I r 

m=0 


cos(m0) + 


X A £pm*r 2 "Tsin(m9) 

1=0 \n=0 J 


(A2) 


while the expressions for multi-valued quantities (e.g., u, 
and u^) assume the following form: 


1 ~ 


M(r,0) = — £A 0l r 2 ' + Ir -f XA ma r 2 " I- 
r " = ) m=l Vn = 0 ) 

cos(m9)+ £ r m ~'( XB ma r 2 * lsin(m0) 
m=l \n-0 J 


(A3) 


As any scalar or Cartesian velocity component is uniquely 
defined at the origin, one can write: 



= 0 


(A4) 


This relation holds, in particular, for u* = u r sin(0)+u ft cos(0), 
where (y,z) plane is oriented perpendicular to the jet axis. 
By taking the derivatives with respect to 0, and requiring 
that the relation holds for any 0, one obtains: 

9u r , 3u ft 

— = u. a„d — „ „0 (AS) 


There is another important constraint on the general form 
of the series expansions for u r and u 0 . If A-j r \ B» r> , Ay 0 * 

and B§ } are the coefficients of the series expansions for u r 


and Ue in (A3), the following relation holds for all i > l : 

A io’ = B !o and B-®' = - A jo (A6) 

By calculating the derivatives with respect to 0 and r of the 
series expansions given by (A2) and (A3) for all operators 
present in the RHS of the governing equations (Al) aid 
taking the limit r -4 0, a new form of the governing 
equations that is valid at the polar axis is obtained. These 
are exact, provided that we can calculate the coefficients 
A mn. Bn». a mo , for all terms present in the RHS of 
(Al). However, for a system of PDE’s with second-oitfer 
radial derivatives, 3s is the case for the Navier-Stokes 
equations, it is sufficient to calculate at most the 
coefficients whose indices m and n vary between 0 and 2. 
For instance, the dilatation operator, which is a scalar 
quantity and should contain only the m=0 mode, can be 
expressed as: 



where the streamwise derivative can be calculated with the 
same method used for points situated away from the polar 
axis. For a complete description of the form of the RHS 
expressions corresponding to the fully compressible flow 
equations including the LES terms, the reader is referred to 
Constantinescu and Lele 17 . 

The last step is to describe how the asymptotic series 
coefficients that are needed to evaluate the RHS in (A 1) are 
computed. All that is required to calculate these 
coefficients accurately is to estimate numerically the first 
and second older radial derivatives of all the variables in 
RHS with the same order of accuracy as for points away 
from the polar axis. To do this, the following algorithm is 
adopted. The computational domain is mapped at every 
x=constant, such that there is no need to specify numerical 
boundary conditions at r=0. The mapping function 
(r, 0) — > (?, 0) is: 


r =r, 

for 

0< 9 < n 

and 

II 

<x> 

0 < r < R 


r f = -r, 

for 

n < 0 < 2n 


0 = 0 - ji 

0< r < R 



The radial derivatives are now taken from -R to R with n=0 
being a regular interior point instead of a ‘numerical* 
boundary point. This is similar to the method proposed by 
Mohseni and Colonius 26 , but in their case points in the 
radial direction were distributed starting with r=Ar/2, instead 
of r=0. 

Once the values of the variables and of their radial 
derivatives are known for all N e directions at r=0 (N 0 is the 
number of points in the azimuthal direction), the relations 
that define the limit of the series expansions at r=0 for the 
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variables and their radial derivatives are used to calculate the 
coefficients in these expansions. For instance, the 
coefficients A\q andB^ that appear in the expansion of 
the radial velocity u,: 

u r = A{q cos(0) + Bjg sin(0) (A9) 

can be calculated by solving the system of equations that is 
obtained by writing (A9) at 9 and 0+rc/2. If N 0 is divisible 
by 8, 11,(9) and u r (0+rc/2) corresponds to modes situated 
apart. As these values are known, the calculation of 
these coefficients is reduced to solving a system of two 
linear equations with two unknowns. To eliminate the bias 
toward a certain direction, one can solve the above system 
for every 0=(2 ti/ N e ) (n-1) with n=l to N 0 and average the 
results to get final values for A\q and bJq. In a similar 
way, the coefficients involved in the expressions for the 
First (and second order derivatives) of u r and u* require 
solving two (and three) systems of two linear equations. 
These coefficients serve also to estimate the limit of terms 
involving the azimuthal derivatives at the polar axis. The 
coefficients in the expansions are unique so we do not have 
to choose any arbitrary directions as was the case when the 
equations were solved in Cartesian coordinates at the 
singularity axis. 
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